*------------------------------------------------------------------------------*
* This program analyzes simulated data from the structural model
*------------------------------------------------------------------------------*

global indata ".../data"

*Import and save annuty factors
import delimited using "$indata/sim_annfactor.csv", clear
rename v1 age
rename v2 annfactor
save "$indata/sim_annfactor", replace

*Import simulated data
import delimited using "$indata/simulated_profiles_s58.csv", clear

*Labels 
rename v1 persid
rename v2 g
rename v3 emp
rename v4 A
rename v5 dca
rename v6 y
rename v7 B_db
rename v8 B_ss
rename v9 mortp
rename v10 trt_flag

*Current age
gen sortord = _n
sort persid sortord trt_flag
bysort persid trt_flag: gen n = _n
gen t = n-6
gen age = 58+t
drop n t

*Merge with ann factors
merge m:1 age using "$indata/sim_annfactor", keepus(annfactor)
keep if _merge==3
drop _merge

*DBwealth
gen dba = B_db*annfactor

*Psuedo persid
egen persid2 = group(persid trt_flag)
xtset persid2 age

*pension accruals
*DC accruals
gen dca_p1 = f.dca
gen dcaccrual = (dca_p1/(1.029)) - dca
*Interpolation error sometimes causes accrual to be<0 
replace dcaccrual = 0 if emp==1 & dcaccrual<0
*DB accruals
gen dba_1 = l.dba
gen dbaccrual = dba - dba_1
gen dbaccrual_f = f.dbaccrual

*deferral fraction
gen defrate = dcaccrual/y
gen mw = (defrate-.011)/1.39 if defrate<.10135
replace mw = (defrate - .03635) if defrate>=.10135
replace mw = 0 if emp==0
replace mw = max(0,mw)

*Deferral amount
gen ern_after_dc = (1-mw)*y


*File to pass to TAXSIM to compute FIITAX liability
preserve

	keep sortord age ern_after_dc 
	
	*file structure for taxsim

	gen taxsimid = sortord   
	gen year = 2010     //2010
	gen state  = 0     
	gen mstat = 1       //1 = single 
	gen page = .        //age
	gen sage = 0       
	gen depx = 0
	gen dep13 = 0
	gen dep17 = 0
	gen dep18 = 0
	gen pwages = .      //annual labor earnings
	gen swages = 0
	gen dividends = 0
	gen intrec = 0
	gen stcg = 0 
	gen ltcg = 0
	gen otherprop = 0
	gen nonprop = 0     
	gen pensions = 0    //pension distributions
	gen gssi = 0        //SS benefits
	gen ui = 0
	gen transfers = 0
	gen rentpaid = 0
	gen proptax = 0
	gen otheritem = 0
	gen childcare = 0
	gen mortgage = 0
	
	replace pwages = ern_after_dc
	replace page = age
	drop age ern_after_dc
	
	*Run the tax calculator
	cd "$indata"
	taxsim35, replace

	*Keep relevant vars
	keep fiitax sortord
	save "$indata/siminput_s58_fiitax", replace
restore

*File to pass to TAXSIM to compute FICA liability (DC contribution taxable)
preserve

	keep sortord age y 
	
	*file structure for taxsim

	gen taxsimid = sortord   
	gen year = 2010     //2010
	gen state  = 0     
	gen mstat = 1       //1 = single 
	gen page = .        //age
	gen sage = 0       
	gen depx = 0
	gen dep13 = 0
	gen dep17 = 0
	gen dep18 = 0
	gen pwages = .      //annual labor earnings
	gen swages = 0
	gen dividends = 0
	gen intrec = 0
	gen stcg = 0 
	gen ltcg = 0
	gen otherprop = 0
	gen nonprop = 0     
	gen pensions = 0    //pension distributions
	gen gssi = 0        //SS benefits
	gen ui = 0
	gen transfers = 0
	gen rentpaid = 0
	gen proptax = 0
	gen otheritem = 0
	gen childcare = 0
	gen mortgage = 0
	
	replace pwages = y
	replace page = age
	drop age y
	
	*Run the tax calculator
	cd "$indata"
	taxsim35, replace

	*Keep relevant vars
	keep fica sortord
	save "$indata/siminput_s58_ficatax", replace
restore


*Merge tax data back 
merge 1:1 sortord using "$indata/siminput_s58_fiitax", keepus(fiitax)
drop _merge
merge 1:1 sortord using "$indata/siminput_s58_ficatax", keepus(fica)
drop _merge

*Total tax liability
gen taxliab = fiitax+fica/2
*replace taxliab = 0 if emp==0

*Total compensation (post-tax)
gen totcomp = (ern_after_dc - taxliab) + dcaccrual + dbaccrual_f if emp==1 
replace totcomp = 0 if emp==0  

*Total compensation defined in retirement if control group case is still employed
sort persid2 age trt_flag
gen totcomp_2 = totcomp
bysort persid age: egen temp = max(totcomp)
replace totcomp_2 = temp if totcomp_2==0
drop temp

*Retirement age
gen temp = age*emp
bysort persid2: egen retage = max(temp) //for each person
*Ret age is the first age out of LF
replace retage=retage+1
gen retage_control = retage*(1-trt_flag) //for the control version
bysort persid: egen retage_all = max(retage_control) //apply the control version to treated group
drop temp retage_control
*Redefine ret age as last working age
replace retage_all = retage_all-1
replace retage = retage-1

*-------------------------------------------------------------------------------
*flow and pdv loss by relative period
*-------------------------------------------------------------------------------
gen k = age-58
gen pdvmtcomp = .
*define counterfactual retage in terms of k (i.e. ret age in control state)
gen kret_all = k if age == retage_all
bysort persid2: egen temp = max(kret_all)
replace kret_all = temp if kret_all==.
drop temp

*Flag individuals working prior to each k
forvalues fk = 0/12{
	gen flag= (emp==1 & k==`fk'-1)
	bysort persid2: egen work_`fk' = max(flag)
	drop flag
}

*Compute elements of the elasticities

matrix eta = J(1,4,.) 

*Frisch

*Employment change
areg emp trt_flag if k==0 & work_0==1, a(persid) cl(persid) 
matrix b = e(b)
matrix V = e(V)

*Numerator of eta and se
scalar temp = b[1,1]/b[1,2]
matrix eta[1,1] = temp
scalar temp =sqrt((1/b[1,2]^2)*V[1,1]+(b[1,1]^2/b[1,2]^4)*V[2,2])
matrix eta[1,2] = temp

*Current compensation change
areg totcomp trt_flag if k==0 & work_0==1, a(persid) cl(persid)
matrix b = e(b)
matrix V = e(V)

*Denominator of eta and se4
scalar temp = b[1,1]/b[1,2]
matrix eta[1,3] = temp
scalar temp =sqrt((1/b[1,2]^2)*V[1,1]+(b[1,1]^2/b[1,2]^4)*V[2,2])
matrix eta[1,4] = temp	

*Output to calculate etas and se's
clear
svmat eta
rename eta1 demp
rename eta2 dempse
rename eta3 dchi
rename eta4 dchise
gen eta = demp/dchi
gen etase = sqrt((1/dchi^2)*dempse^2 + (demp^2/dchi^4)*dchise^2)
export delimited ".../model-output/eta_58.csv", replace





